//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#ifndef __coilgen_hpp__
#define __coilgen_hpp__

#include "ephi.hpp"

#define COILGEN_HEADER_SIZE 128
#define COILGEN_MAGIC "cGe1"
#define COILGEN_BUFFER 4
#define COILGEN_QMOD 1e-6

class CoilGen
{
public:
  CoilGen ();

  prec_t radius;
  prec_t wr;
  prec_t dbl;
  enum { ROUND, SQUARE } shape;
  size_t memory_usage;
  size_t coil_matrix_n;
  bool do_efield;
  bool do_cubic;
  size_t force_grid;

  prec_t geom_radial;
  prec_t geom_axial;

  void set_space (prec_t space)
  {
    dbl = 2 * (radius + (wr * 2 + space)/PREC_SQRT2);
  }

  void generate (const char *path, bool print_info = false, bool show_progress = false);
};

struct coilgen_header
{
  char magic[4];
  unsigned int flags;
  enum eflags {
    IS_HALF                   = 0x00000001,
    IS_COILPAIR               = 0x00000002,
    IS_OPPOSING_COILPAIR      = 0x00000004,
    IS_MIRROR                 = 0x00000008
  };

  unsigned int coilprec_size;
  unsigned int _filler; // to get 8 aligned

  double coil_radius;
  double coil_wr;
  double coilpair_distance;

  unsigned int radial_samples;
  unsigned int axial_samples;
  double geom_radial;
  double geom_axial;
  double geom_axial_add;
};

struct coilgen_blob_header
{
  enum eblobtype {
    IS_BLOB_BMAP = 1,
    IS_BLOB_BMAPEMAP,
    IS_BLOB_CONTACTMAP,
    IS_BLOB_CUBIC_BMAP,
    IS_BLOB_CUBIC_BMAPEMAP,
  };
  unsigned int blob_type;
  unsigned int blob_len_high;
  unsigned int blob_len_low;

  void set_len (size_t len)
  {
    if ( sizeof(size_t) == 4 )
      {
        blob_len_high = 0;
        blob_len_low = len;
      }
#ifndef LONG_IS_32BIT
    else
      {
        blob_len_low = len & 0xffffffff;
        blob_len_high = (len >> 32) & 0xffffffff;
      }
#endif
  }
  bool can_read () const { return sizeof(size_t) >= 8 || !blob_len_high; }
  size_t get_len () const
  {
    return (size_t)blob_len_low
#ifndef LONG_IS_32BIT
        | ((size_t)blob_len_high << 32)
#endif
    ;
  }
};

// float might seem fractionally better on its own, but double is more smooth
// and up to order of magnitude better in relative errors when combining coils
// that is when comparing grids of equal size (MB)
typedef double coilprec_t;
//typedef float coilprec_t;
#if !defined(PREC_FLOAT) && !defined(PREC_DOUBLE)
  #define prec2coilprec(v) prec2double(v)
#else
  #define prec2coilprec(v) (v)
#endif

struct coilgen_bmap
{
  coilprec_t rb;
  coilprec_t ab;
};

struct coilgen_bmapemap
{
  coilprec_t rb;
  coilprec_t ab;
  coilprec_t re;
  coilprec_t ae;
};

struct coilgen_cubic_bmap
{
  typedef coilprec_t matrix_t[4][4];
  matrix_t rb;
  matrix_t ab;
};

struct coilgen_cubic_bmapemap
{
  typedef coilprec_t matrix_t[4][4];
  matrix_t rb;
  matrix_t ab;
  matrix_t re;
  matrix_t ae;
};

class CoilData
{
public:
  CoilData ();
  ~CoilData ();

  void cleanup();

  bool load (const char *path);

  coilgen_bmap *bmap;
  coilgen_bmapemap *bmapemap;
  unsigned char *cmap;
  coilgen_cubic_bmap *cubic_bmap;
  coilgen_cubic_bmapemap *cubic_bmapemap;

  unsigned int flags;

  prec_t coil_radius;
  prec_t coil_wr;

  prec_t coilpair_distance;

  size_t radial_samples;
  size_t axial_samples;
  prec_t geom_radial;
  prec_t geom_axial;
  prec_t geom_axial_add;

  bool has_axial_add;
  bool mirror;

  prec_t cr, ca;
  prec_t rec_cr, rec_ca;
};

#include "polywell_conf.hpp"
#include "ephi.hpp"

template <class IC>
class CoilDataFactory : public CoilFactory
{
public:
  CoilDataFactory (const CoilData& cd) : cd(cd) {}
  virtual ~CoilDataFactory() { }

  bool make (Statics& statics, const vect3d& pos, const vect3d& normal,
      prec_t radius, prec_t wr, prec_t I, prec_t q_per_m) const;

  bool make_coilpair (Statics& statics, const vect3d& pos, const vect3d& normal,
      prec_t distance, prec_t radius, prec_t wr, prec_t I, prec_t q_per_m) const;

private:
  const CoilData& cd;
};

template <class IC>
bool CoilDataFactory<IC>::make (Statics& statics, const vect3d& pos,
    const vect3d& normal, prec_t radius, prec_t wr, prec_t I, prec_t q_per_m) const
{
  if ( cd.coilpair_distance != 0 )
    {
      Ephi::debug (Ephi::DEBUG_CRITICAL, "CoilDataFactory::make: coildata is for a coilpair\n");
      return false;
    }

  if ( cd.coil_radius != radius || cd.coil_wr != wr )
    {
      Ephi::debug (Ephi::DEBUG_CRITICAL, "CoilDataFactory::make: coil radiuses do not match\n");
      return false;
    }

  statics.addStaticElement (new IC(cd, pos, normal, I, q_per_m));
  return true;
}

template <class IC>
bool CoilDataFactory<IC>::make_coilpair (Statics& statics, const vect3d& pos, const vect3d& normal,
    prec_t distance, prec_t radius, prec_t wr, prec_t I, prec_t q_per_m) const
{
  if ( cd.coilpair_distance == 0 )
      return CoilFactory::make_coilpair (statics, pos, normal, distance, radius, wr, I, q_per_m);

  if ( cd.coil_radius != radius || cd.coil_wr != wr )
    {
      Ephi::debug (Ephi::DEBUG_CRITICAL, "CoilDataFactory::make_coilpair: coil radiuses do not match\n");
      return false;
    }

  if ( !cd.mirror )
    {
      Ephi::debug (Ephi::DEBUG_CRITICAL, "CoilDataFactory::make_coilpair: not an opposing coilpair dataset\n");
      return false;
    }

  if ( fabs(cd.coilpair_distance - distance)/distance > 0.000001 )
    {
      Ephi::debug (Ephi::DEBUG_CRITICAL, "CoilDataFactory::make_coilpair: coil distance not matching: %g vs %g\n",
          prec2double(cd.coilpair_distance), prec2double(distance));
      return false;
    }

  statics.addStaticElement (new IC(cd, pos, normal, I, q_per_m));
  Ephi::debug (Ephi::DEBUG_INFO, "CoilDataFactory::make_coilpair: coilpair created\n");
  return true;
}

#endif // __coilgen_hpp__

